
$NOFLOATCALLS
$NODEBUG
$STORAGE:2
c*************************************************************************

	subroutine const(ist,imeth,idef)

c****	subroutine supplies constants used in the runge-kutta
c****	method,in the format required by RKP.

	implicit double precision (a-h,p)
	dimension al(12),b(12,12)
	common /coeffs/al,b,a1,a3,a4,a5,a6,a7,a8,a9,a10,a11,
     &		       a12,a13,er1,er3,er4,er5,er6,er7,er8,
     &		       er9,er10,er11,er12,er13,inum
	common /cons/hmin,hmaxl,hfact,hlim,power

	do 10 j=1,12
	   do 5 i=1,12
	      b(i,j)=0.d0
5	   continue
	   al(j)=0.d0
10     continue

       if(imeth.EQ.1) then
c****  Fourth order runge-kutta method specified by
c****  Felberg, E., COMPUTING,	6(1970)p61-71
c      write(*,*) 'Fourth order Runge-Kutta-Fehlberg method'
       al(1)=1.d0/4.d0
       al(2)=3.d0/8.d0
       al(3)=12.d0/13.d0
       al(4)=1.d0
       al(5)=1.d0/2.d0

       b(1,1)=1.d0/4.d0
       b(1,2)=3.d0/32.d0
       b(2,2)=9.d0/32.d0
       b(1,3)=1932.d0/2197.d0
       b(2,3)=-7200.d0/2197.d0
       b(3,3)=7296.d0/2197.d0
       b(1,4)=439.d0/216.d0
       b(2,4)=-8.d0
       b(3,4)=3680.d0/513.d0
       b(4,4)=-845.d0/4104.d0
       b(1,5)=-8.d0/27.d0
       b(2,5)=2.d0
       b(3,5)=-3544.d0/2565.d0
       b(4,5)=1859.d0/4104.d0
       b(5,5)=-11.d0/40.d0

       er1=1.d0/360.d0
       er3=-384.d0/12825.d0
       er4=-41743.d0/1429560.d0
       er5=1.d0/50.d0
       er6=2.d0/55.d0

       a1=16.d0/135.d0
       a2=0.d0
       a3=6656.d0/12825.d0
       a4=28561.d0/56430.d0
       a5=-9.d0/50.d0
       a6=2.d0/55.d0

       inum=6
       ist=4
       if(idef.eq.0) hmaxl=0.5d0
       hlim1=3.d0
       end if

       if(imeth.EQ.2) then
c****  Fifth order runge-kutta method specified by
c****  Verner J.H., SIAM J. Numer. Anal. V15,(1978),p772.(table 5)
c      write(*,*) 'Fifth order Runge-Kutta-Verner method'
       al(1)=1.d0/18.d0
       al(2)=1.d0/6.d0
       al(3)=2.d0/9.d0
       al(4)=2.d0/3.d0
       al(5)=1.d0
       al(6)=8.d0/9.d0
       al(7)=1.d0

       b(1,1)=1.d0/18.d0
       b(1,2)=-1.d0/12.d0
       b(2,2)=1.d0/4.d0
       b(1,3)=-2.d0/81.d0
       b(2,3)=4.d0/27.d0
       b(3,3)=8.d0/81.d0
       b(1,4)=40.d0/33.d0
       b(2,4)=-4.d0/11.d0
       b(3,4)=-56.d0/11.d0
       b(4,4)=54.d0/11.d0
       b(1,5)=-369.d0/73.d0
       b(2,5)=72.d0/73.d0
       b(3,5)=5380.d0/219.d0
       b(4,5)=-12285.d0/584.d0
       b(5,5)=2695.d0/1752.d0
       b(1,6)=-8716.d0/891.d0
       b(2,6)=656.d0/297.d0
       b(3,6)=39520.d0/891.d0
       b(4,6)=-416.d0/11.d0
       b(5,6)=52.d0/27.d0
       b(1,7)=3015.d0/256.d0
       b(2,7)=-9.d0/4.d0
       b(3,7)=-4219.d0/78.d0
       b(4,7)=5985.d0/128.d0
       b(5,7)=-539.d0/384.d0
       b(7,7)=693.d0/3328.d0

       er1=33.d0/640.d0
       er3=-132.d0/325.d0
       er4=891.d0/2240.d0
       er5=-33.d0/320.d0
       er6=-73.d0/700.d0
       er7=891.d0/8320.d0
       er8=2.d0/35.d0

       a1=57.d0/640.d0
       a3=-16.d0/65.d0
       a4=1377.d0/2240.d0
       a5=121.d0/320.d0
       a7=891.d0/8320.d0
       a8=2.d0/35.d0

       inum=8
       ist=5
       if(idef.eq.0) hmaxl=1.0d0
       hlim1=4.d0
       end if

       if(imeth.EQ.3) then
c****  Seventh order runge-kutta method specified by
c****  Verner J.H., SIAM J. Numer. Anal. V15,(1978),p772.(table 7)
c      write(*,*) 'Seventh order Runge-Kutta-Verner method'
       al(1)=1.d0/4.d0
       al(2)=1.d0/12.d0
       al(3)=1.d0/8.d0
       al(4)=2.d0/5.d0
       al(5)=1.d0/2.d0
       al(6)=6.d0/7.d0
       al(7)=1.d0/7.d0
       al(8)=2.d0/3.d0
       al(9)=2.d0/7.d0
       al(10)=1.d0
       al(11)=1.d0/3.d0
       al(12)=1.d0

       b(1,1)=1.d0/4.d0
       b(1,2)=5.d0/72.d0
       b(2,2)=1.d0/72.d0
       b(1,3)=1.d0/32.d0
       b(3,3)=3.d0/32.d0
       b(1,4)=106.d0/125.d0
       b(3,4)=-408.d0/125.d0
       b(4,4)=352.d0/125.d0
       b(1,5)=1.d0/48.d0
       b(4,5)=8.d0/33.d0
       b(5,5)=125.d0/528.d0
       b(1,6)=-1263.d0/2401.d0
       b(4,6)=39936d0/26411d0
       b(5,6)=-64125.d0/26411d0
       b(6,6)=5520.d0/2401.d0
       b(1,7)=37.d0/392.d0
       b(5,7)=1625.d0/9408.d0
       b(6,7)=-2.d0/15.d0
       b(7,7)=61.d0/6720.d0
       b(1,8)=17176.d0/25515.d0
       b(4,8)=-47104.d0/25515.d0
       b(5,8)=1325.d0/504.d0
       b(6,8)=-41792.d0/25515.d0
       b(7,8)=20237.d0/145800.d0
       b(8,8)=4312.d0/6075.d0
       b(1,9)=-23834.d0/180075.d0
       b(4,9)=-77824.d0/1980825.d0
       b(5,9)=-636635.d0/633864.d0
       b(6,9)=254048.d0/300125.d0
       b(7,9)=-183.d0/7000.d0
       b(8,9)=8.d0/11.d0
       b(9,9)=-324.d0/3773.d0
       b(1,10)=12733.d0/7600.d0
       b(4,10)=-20032.d0/5225.d0
       b(5,10)=456485.d0/80256.d0
       b(6,10)=-42599.d0/7125.d0
       b(7,10)=339227.d0/912000d0
       b(8,10)=-1029.d0/4180.d0
       b(9,10)=1701.d0/1408.d0
       b(10,10)=5145.d0/2432.d0
       b(1,11)=-27061.d0/204120.d0
       b(4,11)=40448.d0/280665.d0
       b(5,11)=-1353775.d0/1197504.d0
       b(6,11)=17662.d0/25515.d0
       b(7,11)=-71687.d0/1166400d0
       b(8,11)=98.d0/225.d0
       b(9,11)=1.d0/16.d0
       b(10,11)=3773.d0/11664.d0
       b(1,12)=11203.d0/8680.d0
       b(4,12)=-38144.d0/11935.d0
       b(5,12)=2354425.d0/458304.d0
       b(6,12)=-84046.d0/16275.d0
       b(7,12)=673309.d0/1636800.d0
       b(8,12)=4704.d0/8525.d0
       b(9,12)=9477.d0/10912.d0
       b(10,12)=-1029.d0/992.d0
       b(12,12)=729.d0/341.d0

       er1=-1.d0/480.d0
       er6=-16.d0/375.d0
       er7=-2401.d0/528000.d0
       er8=2401.d0/132000.d0
       er9=243.d0/14080.d0
       er10=-2401.d0/19200.d0
       er11=-19.d0/450.d0
       er12=243.d0/1760.d0
       er13=31.d0/720.d0

       a1=31.d0/720.d0
       a6=16.d0/75.d0
       a7=16807.d0/79200.d0
       a8=16807.d0/79200.d0
       a9=243.d0/1760.d0
       a12=243.d0/1760.d0
       a13=31.d0/720.d0

       inum=13
       ist=7
       if(idef.eq.0) hmaxl=2.d0
       hlim1=5.d0
       end if

       power=1.d0/dble(ist+1)
       hfact=0.5d0**(1.d0/dble(ist-1))
       hlim=hlim1/hfact
       if(idef.eq.0) hmin=1.d-06
c      write(*,*)' Constants evaluated'
       return
       end
